pacman::p_load(raster, sf, spatstat, tmap)Hands-on Exercise 3: 1st Order Spatial Point Patterns Analysis Methods
1.0 Introduction
1.1 Getting Started
In this hands-on exercise, we will be using the following packages:
maptools for manipulating geographic data
(Note: We will use maptools to convert Spatial objects into ppp format of spatstat)
raster for reading, writing, manipulating, analyzing and modelling of gridded spatial data
(Note: We will use raster to convert image output generated by spatstat into raster format)
sf for handling geospatial data
spatstat for point pattern analysis
(Note: We will use spatstat to perform 1st and 2nd order spatial point patterns analysis and derive kernel density estimation (KDE) layer)
tmap for creating thematic maps such as choropleth and proportional symbol maps,
2.0 Data Acquisition
We will be using 3 datasets in this exercise:
Pre-Schools Location from data.gov.sg
MP14_SUBZONE_WEB_PL from data.gov.sg
MasterPlan2019SubzoneBoundaryNoSeaKML from data.gov.sg
2.1 Extracting Geospatial Data Sets
Following a structure similar to Hands-on Exercise 01, start by creating a new folder labeled Hands-on_Ex03. Within this folder, create a sub-folder named data. Inside the data sub-folder, create one additional sub-folders and rename them geospatial.
Tip: We can shorten MasterPlan2019SubzoneBoundaryNoSeaKML so it is easy to import the dataset.
Let’s rename MasterPlan2019SubzoneBoundaryNoSeaKML to MP19_SUBZONE_BoundaryNoSea.
Unzip MasterPlan2014SubzoneBoundaryWebSHP.zip and place all the unzipped files, PreSchoolsLocation.geojson and MP19_SUBZONE_BoundaryNoSea.kml into the geospatial sub-folder.
3.0 Geospatial Data Handling
3.1 Importing Geospatial Data
In the previous exercises, we have learnt to import geospatial data into RStudio by using st_read() of sf package. Let’s try it now!
childcare_sf <- st_read("data/geospatial/ChildCareServices.geojson")Reading layer `ChildCareServices' from data source
`C:\kt526\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data\geospatial\ChildCareServices.geojson'
using driver `GeoJSON'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84mpsz_sf <- st_read(dsn = "data/geospatial",
layer = "MP14_SUBZONE_WEB_PL")Reading layer `MP14_SUBZONE_WEB_PL' from data source
`C:\kt526\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21sgsz_sf <- st_read(dsn = "data/geospatial/MP19_SUBZONE_BoundaryNoSea.kml")Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source
`C:\kt526\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data\geospatial\MP19_SUBZONE_BoundaryNoSea.kml'
using driver `KML'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY, XYZ
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 843.2 Checking the Contents of A Simple Feature Data Frame
3.2.1 Using st_geometry() to check for inappropriate coordinate systems
st_geometry(childcare_sf)Geometry set for 1925 features
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84
First 5 geometries:st_geometry(mpsz_sf)Geometry set for 323 features
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
First 5 geometries:st_geometry(sgsz_sf)Geometry set for 332 features
Geometry type: MULTIPOLYGON
Dimension: XY, XYZ
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84
First 5 geometries:You might have observed variations in the coordinate systems among the data frames. Recall in Hands-on Exercise 01, it is a common practice to transform original data from geographical coordinate system to projected coordinate system.
childcare_sfandsgsz_sfuses WGS 84 (geographic coordinate system)mpsz_sfuses uses SVY21 (projected coordinate system)
Let’s perform projection transformation on childcare_sf and sgsz_sf using st_transform() of sf package.
childcare3414 <- st_transform(childcare_sf,
crs = 3414)sgsz3414 <- st_transform(sgsz_sf,
crs = 3414)Now, we will display the contents of childcare3414 and sgsz3414.
st_geometry(childcare3414)Geometry set for 1925 features
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 11810.03 ymin: 25596.33 xmax: 45404.24 ymax: 49300.88
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 5 geometries:st_geometry(sgsz3414)Geometry set for 332 features
Geometry type: MULTIPOLYGON
Dimension: XYZ
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 5 geometries:Notice that both childcare3414 and sgsz3414 are now in SVY21 projected coordinate system.
3.2.2 Using st_crs() to check for missing/inaccurate coordinate systems
DIY: Using the appropriate sf function you learned in Hands-on Exercise 2, retrieve the referencing system information of these geospatial data.
To retrieve coordinate reference system from sf object, we need to use st_crs() from sf package.
st_crs(childcare3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]st_crs(mpsz_sf)Coordinate Reference System:
User input: SVY21
wkt:
PROJCRS["SVY21",
BASEGEOGCRS["SVY21[WGS84]",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]st_crs(sgsz3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]Although mpsz_sf is projected in SVY21, the output near the end indicates that EPSG is 9001. This is a wrong EPSG code. The correct EPSG code for SVY21 should be 3414. Let’s change it using st_transform() of sf package.
DIY: Using the method you learned in Lesson 2, assign the correct crs to
mpsz_sf.
mpsz3414 <- st_transform(mpsz_sf, 3414)
st_crs(mpsz3414)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]Notice that the EPSG code is 3414 now.
4.0 Geospatial Visualization
After checking the coordinate reference system of each geospatial data frame, it is also useful for us to plot a map to show their spatial patterns.
DIY: Using the mapping methods you learned in Hands-on Exercise 3, prepare a map as shown below.
tm_shape(sgsz3414) +
tm_polygons() +
tm_shape(mpsz3414) +
tm_polygons() +
tm_shape(childcare3414)+
tm_dots()
Notice that all the geospatial layers are within the same map extend. This shows that their referencing system and coordinate values are referred to similar spatial context. This is very important in any geospatial analysis.
Alternatively, we can also prepare a pin map using tmap_mode() with the view option for interactive mode.
tmap_mode("view")
tm_shape(childcare3414)+
tm_dots()Reminder: Always remember to switch back to plot mode after the interactive map. This is because, each interactive mode will consume a connection. You should also avoid displaying ecessive numbers of interactive maps (i.e. not more than 10) in one RMarkdown document when publish on Netlify.
tmap_mode("plot")5.0 Geospatial Data wrangling
Converting simple feature data frames to sp’s Spatial* classes in geospatial analysis is necessary because many geospatial analysis packages in R, such as “spatial” and “raster,” require input geospatial data in the form of sp’s Spatial* classes.
5.1 Converting sf data frames to sp’s Spatial* class
Let’s convert sf data frames into sp’s Spatial* class using as_Spatial() of sf package.
childcare <- as_Spatial(childcare3414)
mpsz <- as_Spatial(mpsz3414)
sg <- sgsz_sf %>% st_zm() %>% as_Spatial()childcareclass : SpatialPointsDataFrame
features : 1925
extent : 11810.03, 45404.24, 25596.33, 49300.88 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
variables : 2
names : Name, Description
min values : kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>100044</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>44, TELOK BLANGAH DRIVE, #01 - 19/51, SINGAPORE 100044</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td></td> </tr><tr bgcolor=""> <th>NAME</th> <td>PCF SPARKLETOTS PRESCHOOL @ TELOK BLANGAH BLK 44 (CC)</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>349C54F201805938</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20211201093837</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center>
max values : kml_999, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>99982</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>35, ALLANBROOKE ROAD, SINGAPORE 099982</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td></td> </tr><tr bgcolor=""> <th>NAME</th> <td>ISLANDER PRE-SCHOOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>4F63ACF93EFABE7F</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20211201093837</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> mpszclass : SpatialPolygonsDataFrame
features : 323
extent : 2667.538, 56396.44, 15748.72, 50256.33 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
variables : 15
names : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C, REGION_N, REGION_C, INC_CRC, FMEL_UPD_D, X_ADDR, Y_ADDR, SHAPE_Leng, SHAPE_Area
min values : 1, 1, ADMIRALTY, AMSZ01, N, ANG MO KIO, AM, CENTRAL REGION, CR, 00F5E30B5C9B7AD8, 16409, 5092.8949, 19579.069, 871.554887798, 39437.9352703
max values : 323, 17, YUNNAN, YSSZ09, Y, YISHUN, YS, WEST REGION, WR, FFCCF172717C2EAF, 16409, 50424.7923, 49552.7904, 68083.9364708, 69748298.792 sgclass : SpatialPolygonsDataFrame
features : 332
extent : 103.6057, 104.0885, 1.158699, 1.470775 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
variables : 2
names : Name, Description
min values : kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>SUBZONE_NO</th> <td>1</td> </tr><tr bgcolor=""> <th>SUBZONE_N</th> <td>ANAK BUKIT</td> </tr><tr bgcolor="#E3E3F3"> <th>SUBZONE_C</th> <td>BTSZ01</td> </tr><tr bgcolor=""> <th>CA_IND</th> <td>N</td> </tr><tr bgcolor="#E3E3F3"> <th>PLN_AREA_N</th> <td>BUKIT TIMAH</td> </tr><tr bgcolor=""> <th>PLN_AREA_C</th> <td>BT</td> </tr><tr bgcolor="#E3E3F3"> <th>REGION_N</th> <td>CENTRAL REGION</td> </tr><tr bgcolor=""> <th>REGION_C</th> <td>CR</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>4BAD8B2C9CEBF3F2</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20191223152313</td> </tr></table></center>
max values : kml_99, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>SUBZONE_NO</th> <td>9</td> </tr><tr bgcolor=""> <th>SUBZONE_N</th> <td>YISHUN WEST</td> </tr><tr bgcolor="#E3E3F3"> <th>SUBZONE_C</th> <td>YSSZ09</td> </tr><tr bgcolor=""> <th>CA_IND</th> <td>N</td> </tr><tr bgcolor="#E3E3F3"> <th>PLN_AREA_N</th> <td>YISHUN</td> </tr><tr bgcolor=""> <th>PLN_AREA_C</th> <td>YS</td> </tr><tr bgcolor="#E3E3F3"> <th>REGION_N</th> <td>NORTH REGION</td> </tr><tr bgcolor=""> <th>REGION_C</th> <td>NR</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>95C11920195B86C7</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20191223152313</td> </tr></table></center> Notice that the geospatial data have been converted into their respective sp’s Spatial* classes now.
5.2 Converting the Spatial* class into generic sp format
spatstat requires the analytical data in ppp object form. There is no direct way to convert a Spatial* classes into ppp object. We need to convert the Spatial classes* into Spatial object first.
childcare_sp <- as(childcare, "SpatialPoints")
mpsz_sp <- as(mpsz, "SpatialPolygons")
sg_sp <- as(sg, "SpatialPolygons")childcare_spclass : SpatialPoints
features : 1925
extent : 11810.03, 45404.24, 25596.33, 49300.88 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs mpsz_spclass : SpatialPolygons
features : 323
extent : 2667.538, 56396.44, 15748.72, 50256.33 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs sg_spclass : SpatialPolygons
features : 332
extent : 103.6057, 104.0885, 1.158699, 1.470775 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs 5.3 Converting the generic sp format into spatstat’s ppp format
Now, we will use as.ppp() function of spatstat to convert the spatial data into spatstat’s ppp object format.